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ABSTRACT 

The  Smart  Weapons  Test  Range  (SWTR)  lies  within  the  Yuma  Proving  Ground  (YPG),  Arizona.  SWTR 
is  a  new  facility  constructed  specifically  for  the  development  and  testing  of  futuristic  intelligent  battle¬ 
field  sensor  networks.  In  this  paper,  results  are  presented  for  an  extensive  high-resolution  geophysical 
characterization  study  at  the  SWTR  site  along  with  validation  using  3-D  modeling.  In  this  study,  several 
shallow  seismic  methods  and  novel  processing  techniques  were  used  to  generate  a  3-D  grid  of  earth  seis¬ 
mic  properties,  including  compressional  (P)  and  shear  (S)  body-wave  speeds  (Vp  and  Vs),  and  their 
associated  body-wave  attenuation  parameters  (Qp,  and  Qs).  These  experiments  covered  a  volume  of  earth 
measuring  1500  m  by  300  m  by  25  m  deep  (11  million  cubic  meters),  centered  on  the  vehicle  test  track  at 
the  SWTR  site.  The  study  has  resulted  in  detailed  characterizations  of  key  geophysical  properties.  To  our 
knowledge,  results  of  this  kind  have  not  been  previously  achieved,  nor  have  the  innovative  methods 
developed  for  this  effort  been  reported  elsewhere.  In  addition  to  supporting  materiel  developers  with 
important  geophysical  information  at  this  test  range,  the  data  from  this  study  will  be  used  to  validate 
sophisticated  3-D  seismic  signature  models  for  moving  vehicles. 

Each  of  the  four  material-property  volumes  (Vp,  Vs,  Qp,  and  Qs)  was  constructed  by  interpolating  2-D 
seismic  measurements  made  along  survey  lines  optimized  for  the  physical  properties  at  the  SWTR  site. 
The  geostatistical  properties  of  the  data  guided  the  interpolation  of  data  points  between  survey  lines  and 
established  confidence  limits  around  each  interpolated  value.  Ground  truth  was  accomplished  through 
cross-hole  seismic  measurements  and  borehole  logs.  Surface  wave  and  refraction  acquisition  methods 
were  used  to  acquire  most  raw  data  for  this  study.  In  addition  to  standard  reflection  and  refraction  data 
analysis  procedures,  we  also  applied  turning  ray  tomography  and  surface  wave  analysis.  A  variety  of 
seismic  energy  sources  (including  vibroseis,  accelerated  weight  drop,  and  high  frequency  impulses  from 
projectile  and  explosive  shots)  and  recording  array  configurations  were  considered  for  optimizing  the 
quality  of  each  analysis. 

Inversion  of  the  specific  portions  of  the  seismic  wavefield  was  key  to  the  success  of  the  characterization 
effort.  Tomographic  inversion  produced  P-wave  velocity  matrices  with  the  greatest  resolution  and  finest 
sampling  without  the  necessity  for  assumptions  about  layer  properties.  Standard  refraction  analysis  was 
used  to  confirm  the  tomographically  defined  P-wave  velocities  and  to  map  prominent  layers.  Inversion  of 
surface-wave  dispersion  curves  from  multichannel  data  resulted  in  optimal  estimates  of  the  S-wave 
velocity  field.  Prior  to  this  study,  no  numerically  rigorous  methods  had  been  documented  for  estimating  Q 
for  near-surface  materials.  Geophysical  properties  at  the  site  appear  to  be  anisotropic,  reflecting  local 
trends  in  near-surface  geology.  For  example,  spatial  continuity  in  Vs  can  be  envisioned  as  an  ellipsoid 
whose  semiaxes  are  1 19  m,  95  m,  and  9.5  m  in  length.  The  longest  axis,  representing  the  range  of  greatest 
spatial  continuity,  is  oriented  W23°N  and  the  shortest  axis  is  vertical.  Block  kriging  was  used  to  estimate 
averages  of  geophysical  properties  within  0.5  x  0.5  x  0.5  m  cells  and  the  uncertainties  associated  with 
each  block  in  the  volume.  Weights  used  in  kriging  are  a  function  of  distance,  direction,  and  rate  of  spatial 
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change  observed  in  the  data.  The  noted  uncertainties  tended  to  reach  their  largest  values  for  cells  located 
farthest  from  measurement  points,  but  rarely  exceeded  10%  of  the  estimated  values  (Vs,  Vp,  Qp,  and  Qs). 

Forward  modeling  was  conducted  on  the  Yuma  3-D  model  using  a  3-D  8th  order  finite  difference  scheme 
that  includes  topography.  The  model  had  a  dimension  401  x  151  x51  cells  with  0.5  meter  spacing.  The 
time  step  used  was  0.0001  seconds.  Due  to  the  size  of  the  model,  DoD  supercomputers  were  used  in  the 
computation.  A  vertical  point  source  was  used  with  peak  pressure  at  0.015  seconds.  The  elastic  (infinite 
Q)  model  was  computed  for  reference  and  to  bench  mark  the  attenuation.  Results  are  compared  with 
Yuma  field  data.  After  scaling  the  source  magnitude  good  agreement  is  found  in  amplitude,  arrival  times, 
amplitude  decay  rates  and  frequency  spectra  of  dispersed  phases. 

1.  INTRODUCTION 

Development  of  battlefield  sensor  networks  intended  to  classify,  track,  and  monitor  ground  vehicle  traffic 
in  forward  operation  areas  will  become  increasingly  critical  to  the  modem  Army’s  intelligence  mission. 

It  is  expected  that  these  networks  will  rely  heavily  on  seismic  sensing  capabilities  and  must  function  with 
a  high  degree  of  reliability  across  a  wide  range  of  geologic  settings.  Discriminating  and  classifying  vari¬ 
ous  types  of  surface  vehicles  based  on  their  unique  seismic  characteristics  or  signature  requires  a  priori 
knowledge  of  each  target  vehicle’s  characteristics  and  information  about  the  local  near-surface  earth 
properties  that  influence  the  propagation  of  the  seismic  wavefield.  Confidently  estimating  key  near¬ 
surface  earth  material  properties  in  hostile  areas  will  ultimately  require  interfacing  near-surface  geo¬ 
physics  and  remote  sensing  technologies.  Earth  seismic  properties  are  critical  components  to  wavefield 
simulations  and  ultimately  must  be  available  for  incorporation  into  sensors  designed  for  deployment  in 
areas  unavailable  for  occupation. 

Effective  seismic  sensing  technologies  deployed  in  forward  operational  areas  will  depend  on  direct  or 
inferred  site -specific  knowledge  of  the  earth’s  physical  properties.  Wavefield  simulations  and  verification 
of  smart-system  performance  requires  calibration  in  specific  type  localities.  Such  perfomiance  evaluations 
require  accurate  estimates  of  real  seismic  properties  in  three  dimensions,  with  a  level  of  3-D  detail  con¬ 
sistent  with  site -specific  geology  and  computational  constraints.  Accuracy  and  resolution  of  the  metho¬ 
dologies  used  to  extract  critical  earth  seismic  properties  form  the  basis  for  confidence  ratings  and  boun¬ 
dary  conditions  required  to  develop  procedures  for  vehicle -type  identification  from  seismic  wavefield 
analysis.  Correlation  of  seismic  wavefield  models  with  real  seismograms  in  specific  geologic  settings 
provides  both  the  mechanism  to  quantify  the  performance  of  simulations  and  discrimination  algorithms, 
and  the  best  means  to  refine  various  approaches  of  extracting  seismic  properties  from  field  seismic 
wavefield  measurements. 

Resolution  requirements  of  seismic  wavefield  propagation  models  and  real  earth  properties  dictate  the 
size  and  orientation  of  each  grid  cell  within  the  specified  earth  volume.  Characteristics  of  this  volume 
determine  the  propagation  of  the  modeled  seismic  wavefield  used  for  quantitative  comparisons  with  field- 
measured  vehicle  signature  data.  By  incorporating  geostatistical  analyses  into  the  estimation  of  geo¬ 
physical  properties  in  each  cell  of  the  volume,  the  estimations  can  be  optimized  and  the  estimation  pro¬ 
cess  made  more  efficient.  Geostatistics  focuses  on  the  rate  and  direction  of  changes  in  geophysical 
characteristics  of  the  near  surface  and  relates  these  changes  to  the  spatial  arrangement  of  the  sample 
points.  Designing  the  distribution  of  cells  in  the  model,  interpolating  between  sample  points,  and  con¬ 
structing  confidence  levels  for  geophysical  properties  in  the  3-D  volume  is  best  done  using  geostatistical 
procedures. 

Key  seismic  propagation  parameters,  including  seismic  body-wave  velocities  (Vp,  Vs)  and  estimates  of 
seismic  energy  loss  in  terms  of  Qp  and  Qs,  must  be  determined  for  each  cell.  In  addition  to  these  key  prop¬ 
erties,  earth  models  must  also  preserve  local  topographic  variation,  layer  positioning,  and  thicknesses  to 
be  realistic.  Since  seismic  properties  can  be  defined  using  data  acquired  via  a  variety  of  geophysical 


methodologies  and  approaches,  it  is  necessary  to  evaluate  each  to  insure  the  most  direct,  accurate,  expe¬ 
dient,  and  (laterally  and  vertically)  continuous  approach  is  employed  for  each  unique  geologic  setting. 
Evaluations  of  propagation  patterns  and  the  imaging  potential  of  different  seismic  methods  focused  on 
P-and  S-  body  waves  (reflected,  direct,  and  refracted)  as  well  as  multi-modal  surface  waves  (inversion  of 
dispersive  attributes).  Estimates  of  these  seismic  properties  throughout  the  3-D  earth  model  can  effec¬ 
tively  be  bounded  by  geostatistical  analyses. 

Optimal,  and  preferably  redundant,  independent  estimates  of  Vp,  Vs,  Qp,  and  Qs  are  necessary  for  each 
subsurface  cell.  Acquiring  seismic  data  using  parameters  and  equipment  targeting  specific  parts  of  the 
wavefield  permit  the  most  accurate  estimation  of  Vp  and  Vs.  Estimating  Q  compatible  with  the  rigors  of 
near-surface  wavefield  modeling  has  received  only  scant  notice  in  the  published  scientific  literature  (Jeng, 
et  al.,  1999).  Considering  the  fundamental  dependence  of  surface  wave  energy  on  P-  and  S-  wave  charac¬ 
teristics  of  earth  materials,  determining  Rayleigh  wave  attenuation  as  a  function  of  frequency  allows 
estimations  of  Q  for  both  S-  and  P-type  waves.  Extending  frequency-dependent  attenuation  into  depth 
requires  innovative  approaches  to  sampling  and  inversion. 

The  ability  to  accurately  predict  wave  propagation  in  complex  geology  can  have  a  huge  impact  on  both 
sensor  system  development  time  and  expense.  Previously,  the  limiting  factor  in  3-D  seismic  wave  model¬ 
ing  has  been  that  3-D  representations  of  geology  make  the  computational  task  of  forward  modeling  too 
expensive.  Recent  leaps  in  computational  power  of  the  supercomputers  make  this  a  possibility.  While  the 
physics  regarding  wave  propagation  is  well  known,  validation  of  the  computer  algorithms  and  the  earth 
model  are  necessary  to  make  virtual  sensor  system  development  on  computers  a  viable  option.  With 
validated  models  there  is  no  limit  to  model  scenarios  or  sensor  placement,  making  virtual  sensor  system 
development  a  possibility  and  saving  millions  of  dollars  in  system  development. 

Up-stream  support  of  such  a  virtual  sensing  system  necessitated  a  set  of  carefully  designed,  near-surface 
geophysical  experiments  using  state-of-the-art  high-resolution  seismic  exploration  methods  combined 
with  geostatistical  analysis  techniques  to  provide  essential  earth  properties  for  the  desert  setting  of  YPG. 
The  primary  objective  of  these  experiments  was  to  define  3-D,  micro-geophysical  variations  over  a  well- 
defined  portion  of  the  SWTR  site  at  YPG.  Of  particular  concern  was  the  mapping  the  seismic  earth 
properties  that  most  significantly  impact  the  multi-mode  propagation  of  the  seismic  surface  waves. 

2.  SITE  CHARACTERIZATION  USING  GEOPHYSICS 

Historically,  when  geophysics  has  been  used  for  site  characterization,  the  attempt  involves  measurement 
and  incorporation  of  data  from  at  most  two  or  perhaps  three  different  tools  or  techniques,  each  responding 
differently  to  different  earth  properties.  Rendering  of  these  data  into  quantitative  information  about  the 
geophysical  characteristics  of  the  subsurface  has  generally  stopped  once  first  order  correlations  with 
specific  property  changes  have  been  made  with  “reasonable”  confidence.  For  example,  if  shallow  seismic 
is  one  of  the  tools,  interpretations  typically  have  focused  on  qualitative  assessment  of  subsurface  layer 
topography  rather  than  on  quantitative  measurements  (Clement,  et  al.,  1997;  Pullan  and  Hunter,  1990; 
Lankston,  1990).  The  resulting  acoustically  derived  subsurface  structure  and  velocity  maps  are  incorpo¬ 
rated  into  geologic  and  hydrologic  models  constructed  from  borehole  data.  These  simplistic  models  are 
then  used  for  groundwater  monitoring  purposes  in  an  intuitive,  experience-based  manner  (Steeples  and 
Miller,  1 990;  Miller  and  Xia,  1 999). 

Considering  the  wealth  of  information  contained  in  the  seismic  wavefield,  seismic  measurements  or  imag¬ 
ing  data  have  been  underutilized  for  site  characterization  (Steeples,  et  al.,  1995).  Surface  seismic  tech¬ 
niques  have  been  limited  almost  exclusively  to  routine  mapping  and  delineation  of  subsurface  structures, 
layer  topography,  anomalies,  and  stratigraphic  changes  (Jongerius  and  Helbig,  1988;  Miller,  et  al.,  1989; 
Goforth  and  Hayward,  1992;  Miller,  et  al.,  1995;  Shtivelman,  et  al.,  1998;  Guo  and  Liu,  1999;  Stokoe, 


et  al.,  1994;  Michaels,  1999).  In  many  instances,  different  earth  properties  (Vp,  Vs,  Qp,  Qs,  layer  orienta¬ 
tions,  and  thicknesses)  could  be  estimated  from  different  parts  of  the  seismic  wave  type,  providing  the 
potential  for  redundant  measurement  if  several  techniques  have  been  applied  coincidently. 

3.  GEOLOGIC  SETTING 

YPG  is  in  the  Sonoran  Desert  of  the  Basin  and  Range  Province  in  southwestern  Arizona  (Millet  and 
Barnett,  1970).  In  general,  this  area  is  tectonically  characterized  by  crustal  extension  that  has  produced 
sequences  of  grabens  separating  fault-bounded  mountains  (horsts).  These  bedrock  uplifts  (mountains)  are 
separated  by  broad  valleys  filled  with  alluvial  and  colluvial  sediments  eroded  from  those  mountains.  Sedi¬ 
ments  that  fill  the  basins  have  grain  sizes  that  range  from  clay  to  boulders.  The  SWTR  site  is  within  the 
KOFA  firing  range  where  the  ground  surface  is  characterized  by  extensive  outwash  networks  separated 
by  fan  deposits  and  areas  of  desert  pavement.  The  sparse  vegetation  is  mostly  clustered  around  washes. 

The  ground  surface  at  the  SWTR  site  slopes  gently  southward  with  distinct  1  to  2  m  high  ridges  that 
meander  generally  north  to  south  across  the  northern  half  of  the  site.  Sediment  particles  within  the 
0.5  km2  test  area  vary  in  size  from  gravel  to  clay,  with  locations  in  the  site  where  sorting  into  distinct 
surface  patches  is  evident.  Geological  logs  from  borings  near  the  center  of  the  SWTR  site  suggest  that 
the  first  100  m  of  sediment  are  composed  predominantly  of  silt,  gravel,  and  sand  (Jones,  2001).  From  0 
to  7  m  depth,  the  materials  are  silt  with  sand  and  minor  amounts  of  gravel.  Between  7  and  28  m  depth, 
the  sediments  are  mostly  gravel  with  some  suggestions  of  boulders.  From  28  to  over  65  m  depth,  the 
sediments  are  predominantly  silts  and  sands  with  minor  amounts  of  gravel. 

4,  TESTING  VARIOUS  SEISMIC  METHODS 

Establishing  feasibility  and  appropriate  acquisition  design  specifications  for  quantifying  the  3-D  seismic 
properties  of  the  test  area  mandated  a  unique  combina¬ 
tion  of  seismic  walkaway  testing  with  geostatistical 
analyses.  Test  data  were  acquired  using  a  2'/2-D  radial- 
fan  spread  configuration,  which  from  a  geostatistical 
perspective  best  provided  estimates  of  geophysical 
directionality  in  the  test  area  (Figure  1).  To  accommo¬ 
date  the  requirements  of  this  novel  testing  geometry 
and  procedure,  constraints  on  depth  of  interest,  logis¬ 
tical  limitations,  and  practicality  issues  were  also  con¬ 
sidered  when  deploying  the  radial  source  and  receiver 
array.  Data  acquired  along  each  of  three  radial  pro¬ 
files  were  used  to  estimate  the  orientation  and  extent 
of  spatial  continuity  of  geophysical  properties.  The 
radial  pattern  was  aligned  with  and  near  the  center  of 
the  test  track  at  the  SWTR  site  and  in  close  proximity 
to  the  test  borings  (Figure  2).  These  walkaway  noise 
tests  were  the  basis  for  determining  the  accuracy  and 
effectiveness  of  each  seismic  method,  as  well  as 
estimating  the  optimal  orientation  and  spacing  of  the 
seismic  array. 

.  .  Figure  1.  Radial  walkaway  spread  deployed  to  in ves- 

A  variety  of  sources  and  receivers  was  evaluated  using  tigate  the  antisotropy  of  geophysical  properties  at  the 

this  radial  spread  pattern.  Receiver  separation  for  the  SWTR  site  From  a  geostatistical  perspective  this 
240  receivers  making  up  a  line  in  the  fan  pattern  was  source/receiver  configuration  best  provided  esti- 

0.5  m  or  1 .0  m  (depending  of  the  energy  mode  being  mates  of  geophysical  directionality  in  the  test  area. 


Figure  2.  SWTR  site  relative  to  the  state  of  Arizona. 
Profile  lines  as  defined  by  the  high-resolution  GPS 
survey.  Relative  grid  orientations,  testing  areas,  and 
well  locations  are  identified. 


analyzed)  with  source  stations  located  240  m,  120  m, 
and  1  m  off  each  end  and  centered  on  the  spread. 

With  this  configuration  all  critical  offsets  and  propa¬ 
gation  directions  could  be  fully  sampled  for  each 
source  configuration  and  energy  type.  Four  different 
receivers  were  used,  depending  on  the  portion  of  the 
wavefteld  being  targeted  and  the  dominant  mode  of 
seismic  waves  being  generated  (14  Hz  horizontal, 

4.5  Hz  vertical,  10  Hz  vertical,  and  40  Hz  vertical). 

A  total  of  nine  different  sources  was  tested  at  appro¬ 
priate  offsets.  These  include  20  lb  and  12  lb  sledge¬ 
hammers,  12-  and  8-gauge  downhole  shotguns, 

30.06  and  50-cal  downhole  guns,  LSS6  land  airgun, 
Rubberband  Assisted  Weight  Drop  (RAWD),  and  1V1 
minivib. 


Analysis  of  these  data  for  high-resolution  P-wave  or 
S-wave  seismic  reflections  using  both  conventional 
and  high-resolution  methods  (Steeples  and  Miller, 
1990),  as  well  as  adaptations  made  specifically  for 
these  data  did  not  reveal  events  that  could  be  confi¬ 
dently  and  consistently  inteipreted  as  reflections. 
Model  reflection  hyperbolae  generated  from  borehole 
data  were  used  as  templates  for  evaluating  coherent 
arrival  patterns.  In  light  of  the  extensive  testing  and 
analysis  of  these  data,  it  is  unlikely  that  reflected 
energy  having  significant  resolution  potential  has 
been  produced,  much  less  recorded,  at  this  site.  The 
focus  of  the  reflection  data  acquisition  and  analysis 
was  the  generation  and  identification  of  the  high- 
resolution  reflections  (>100  Hz  P-wave  and  60  Hz 
S-wave)  necessary  for  successful  use  of  this  method. 


Surface  waves  traditionally  have  been  viewed  as 
noise  in  multichannel  seismic  data  collected  to  image 
targets  for  shallow  engineering,  environmental,  and 
groundwater  purposes  (Steeples  and  Miller,  1990). 
Recent  advances  in  the  use  of  surface  waves  for  near¬ 
surface  imaging  combine  spectral  analysis  techniques 
(SASW),  developed  for  civil  engineering  applications 
(Nazarian,  et  al.,  1983),  with  multitrace  reflection 
technologies  developed  for  near  surface  (Schepers, 
1975)  and  petroleum  applications  (Glover,  1959). 

The  combination  of  these  two  uniquely  different 
approaches  to  seismic  imaging  of  the  shallow  sub¬ 
surface  permits  non-invasive  estimation  of  S-wave 
velocities  and  delineation  of  horizontal  and  vertical 


4.1  Reflection 


4.2  Surface  Wave  Inversion 


variations  in  near-surface  material  properties  based  on  changes  in  these  velocities  (MASW)  (Park,  et  al., 
1996;  Xia,  et  ah,  1999;  Park,  et  ah,  1999). 

Extending  this  imaging  technology  to  include  lateral  variations  in  lithology  as  well  as  tunnel  and  frac¬ 
ture  detection,  bedrock  mapping,  and  delineation  of  subsidence  in  karst  terrains  has  required  a  unique 
approach  that  incorporates  SASW,  MASW,  and  CMP  methods.  By  integrating  these  techniques,  2-D 
S-wave  velocity  profiles  of  the  subsurface  can  be  generated.  Estimating  the  dispersion  curve  from  up  to 
60  closely  spaced  receiving  channels  calculated  every  1  to  5  m  along  the  ground  surface  enhances  the 
signal  and  results  in  a  unique,  relatively  continuous  view  of  shallow  subsurface  S-wave  velocity  charac¬ 
teristics.  This  highly  redundant  method  consistently  produces  S-wave  velocity  estimates  that  are  within 
15%  of  borehole  measurements  (Xia,  et  ah,  2000).  Continuous  oversampling  (redundant  sampling)  of  the 
subsurface  minimizes  the  likelihood  that  irregularities  resulting  from  erratic  dispersion  curves  or 
inversion  will  corrupt  the  analysis. 

Estimates  of  Q  flow  naturally  from  analyses  of  surface  wave  attenuation,  specifically  from  Rayleigh  wave 
attenuation  (Xia,  et  ah,  2001).  Q  values  for  the  shallow  subsurface  were  estimated  using  three  related  but 
different  approaches  ('AX  from  velocity  and  ar(f),  single-value  average  of  near-surface  layer  based  on  'AX 
approximation,  and  inversion  of  ar(f)  using  Vs  derived  from  MASW),  each  solving  for  Qs  and  estimating 
Qp  based  on  Qp~  2QS.  Inverting  a,(f)  to  solve  for  Qs  using  Vs  derived  from  MASW  is  significantly  more 
stable  and  accurate  (reduces  error  [>25%])  than  commonly  employed  AVO  methods  (Jin,  et  ah,  2000). 


Surface-wave  data  from  the  radial  test  spread  were  of  sufficient  bandwidth  and  signal-to-noise  ratio  to 
allow  consistent  and  convergent  inversion  of  dispersion  curves  to  S-wave  velocities.  Analysis  of  the  test 
data  allowed  selection  of  sources,  receivers,  and  spread/offset  ranges  optimal  for  sitewide  investigations. 
Since  source,  receiver,  and  spectral  requirements  for  optimal  surface  wave  data  were  different  from  those 


required  to  collect  optimal  near-surface  reflec¬ 
tion  data,  if  both  reflection  and  surface  wave 
methods  had  provided  useful  information  it 
would  have  been  necessary  to  acquire  at  least 
two  unique  data  sets  along  each  profile. 

4.3  Refraction/Tomography 

Direct  and  refracted  P-wave  and  S-wave 
arrivals  were  analyzed  using  conventional 
methods  (Palmer,  1981;  Haeni,  1986;  Lankston, 
1990)  and  inversion  techniques  (Scott,  1977; 
Schneider,  et  ah,  1992;  Ivanov,  et  ah,  2000). 
Use  of  direct  and  refracted  arrivals  for  mapping 
distinct  velocity  contrasts  between  layers  has 
been  in  routine  use  for  everything  from  crustal 
seismic  research  (Steinhart  and  Meyer,  1 96 1 )  to 
shallow  groundwater  studies  (Haeni,  1978). 
Well-documented  limitations  of  the  technique 
(Soske,  1954;  Sander,  1978)  associated  with 
velocity  inversions  and  hidden  layers  limited 
refraction  analysis  at  this  site.  A  pronounced 
velocity  inversion  in  both  P-  and  S-wave  bore- 
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hole  measurements  at  30  m  constrains  classic  Figure  3.  Vp  and  Vs  profiles  from  the  boreholes  located  at 


refraction  data  analysis  (Figure  3).  Methods  about  the  mid-way  point  of  the  test  tract  (Figure  2)  (Block, 
exist  for  approximating  solutions  when  physical  2001).  Measurements  based  on  cross-hole  configuration. 


conditions  (such  as  the  velocity  inversion  observed  at  the  SWTR  site)  violate  assumptions  of  the  refrac¬ 
tion  method  (Mooney  1981;  Redpath  1973),  but  lack  the  robust  nature  of  tomographic  techniques  when 
dealing  with  these  situations. 

Tomography  has  been  used  to  solve  many  subsurface  problems  (Peterson  et  al,  1985;  Cottin,  et  al.,  1986; 
Lytle  and  Dines,  1980;  Kilty  and  Lange,  1990).  A  tomographic  technique  (Joint  Analysis  of  Surface 
Waves  and  Refraction  [JASR])  incoiporating  inversion  of  first  arrivals  with  an  initial  model  from  S-wave 
velocity  profiles  from  surface  wave  data  (Ivanov,  et  al.,  2000)  provided  2-D  Vp  sections  consistent  with 
borehole  measurements.  The  simplicity  of  acquisition  and  computations  made  tomographic  analysis 
especially  attractive  for  velocity  estimation  using  data  optimally  acquired  for  surface-wave  analyses. 
Perhaps  equally  important  is  the  inherent  subdivision  of  the  earth  into  cells,  fundamental  to  tomographic 
analysis.  A  true  advantage  to  the  JASR  method  used  here  is  its  ability  to  overcome  the  non-uniqueness 
problem  inherent  in  conventional  refraction  and  refraction  tomography  methods  (Ivanov,  et  al.,  2000). 
This  approach  increases  the  detail  in  resulting  images  and  therefore  improves  the  apparent  resolution,  and 
since  calculations  of  Vp  are  based  on  a  cellular  approach,  correlations  with  the  Vs  grid  cells  calculated  by 
MASW  are  straightforward. 

First  arrival  analysis  of  the  radial  pattern  walkaway  test  data  was  limited  to  the  upper  30  m  and  produced 
a  two-layer  Vp  model  using  conventional  refraction  methods  and  a  10-layer  Vp  cell  model  using  the  JASR 
method.  The  continuity  in  transitions  between  layers  inherent  in  tomographic  analysis  provides  a  more 
meaningful  cross-section  for  incorporation  into  grid  cells  containing  other  seismic  properties.  Both  Vp 
models  were  produced  for  the  sitewide  2'A-D  grid,  but  only  the  JASR  data  were  incorporated  into  the 
wavefield  simulations.  Secondary  analysis  of  the  radial  pattern  walkaway  data  indicated  that  parameters 
and  equipment  used  to  optimally  acquire  surface  wave  data  for  the  SWTR  site  were  well  within  the  toler¬ 
ances  for  both  conventional  refraction  and  JASP  first  arrival  analysis.  Using  the  same  data  set  for  several 
analyses  minimized  the  field  data  collection  effort  and  enhanced  the  consistency  in  estimates  of  different 
seismic  properties  for  each  cell. 

5.  GEOSTATISTICAL  ANALYSIS  AND  GRID  DESIGN 

A  major  objective  of  the  walkaway  test  was  to  determine  if  there  were  significant  anisotropies  in  geo¬ 
physical  properties  at  the  test  site,  possibly  as  a  result  of  local  trends  in  the  geology.  The  radial  design  of 
the  seismic  lines  provided  measurements  of  geophysical  properties  along  three  different  orientations.  By 
computing  the  semivariance  along  these  three  lines,  a  directional  estimate  of  the  range  of  a  geophysical 
property  such  as  Vp  could  be  made.  (The  range  defines  the  spatial  limit  of  partial  dependence  between 
observations;  beyond  the  range,  observations  are  statistically  independent  and  values  at  one  location  are 
unrelated  to  values  at  other  locations  which  lie  beyond  the  range.)  In  three  dimensions,  the  range  can  be 
envisioned  as  an  ellipsoidal  envelope.  The  three  walkaway  lines  can  be  regarded  as  vectors  extending 
from  the  center  of  the  ellipsoid  to  the  margins  of  the  envelope.  The  direction  cosines  of  the  vectors 
express  the  orientations  of  the  walkaway  lines  and  their  lengths  are  proportional  to  the  ranges  of  the 
geophysical  properties  along  the  lines.  These  relationships  can  be  expressed  in  a  3  x  3  matrix  whose 
eigenvalues  and  eigenvectors  represent  the  major  and  minor  axes  of  the  envelope,  and  which  define  the 
direction  and  extent  of  anisotropy.  At  the  SWTR  site,  the  direction  having  the  greatest  range  is  oriented 
N23°W  with  a  range  of  1 1 8.9  m.  The  range  in  the  perpendicular  direction  is  95. 1  m.  The  range  in  the 
vertical  direction  is  9.5  m.  These  properties  allow  the  design  of  an  optimal  sampling  grid  for  the  main 
phase  of  data  collection. 

The  sampling  plan  used  at  the  SWTR  site  consists  of  parallel  lines  oriented  in  the  direction  of  most  rapid 
change  in  geophysical  properties.  Observations  collected  along  these  lines  have  the  greatest  sampling 
density  (1.2  m)  in  the  direction  of  most  rapid  change  in  seismic  properties.  The  spacing  between  lines 
(50  m)  is  much  greater,  but  seismic  properties  are  more  persistent  in  this  direction.  The  line  spacing  of 


50  m  is  such  that  there  is  complete  overlap  in  range  from  one  line  to  the  next,  so  every  estimated  location 
within  the  array  is  partially  related  to  at  least  two  seismic  lines.  Either  through  good  fortune  or  foresight 
by  the  designers  of  the  SWTR  site,  the  direction  of  maximum  range  in  geophysical  properties  is  approxi¬ 
mately  parallel  to  the  test  track,  which  greatly  eased  the  layout  of  the  seismic  array. 

6.  ACQUISITION  OF  2V2-D  GRID 

Geophysical  and  geo  statistical  analyses  of  walkaway  test  data  using  seismic  velocities  measured  from 
cross-hole  studies  for  control  formed  the  basis  for  selecting  optimal  sources  and  receivers,  source  and 
receiver  spacings,  line  spacings,  and  the  orientation  of  the  major  axis  of  the  2'A-D  rectangular  array.  The 
primary  motivation  for  the  radial  walkaway  test  was  to  allow  examination  of  a  wide  range  of  methods  and 
equipment  to  determine  the  optimum  parameters  and  survey  design,  based  on  geostatistical  and  geo¬ 
physical  evaluations. 

After  geostatistical  analysis  of  the  radial  test  data,  a  214-D  grid  with  the  long  axis  oriented  parallel  to  the 
test  track  was  determined  to  be  optimum  for  local  geologic  conditions  (Figure  2).  Inline  spacing  of 
geophones  was  1.2  m  with  short  axis  lines  300  m  long  and  separated  by  50  m.  Tie  lines  were  separated 
by  100  m,  parallel  to  the  test  track,  and  50  m  offset  from  the  center  of  the  300  m  spreads  and  test  track. 
This  40-spread  configuration  resulted  in  a  uniform  214-D  grid  within  the  1.5  km  by  300  m  study  area. 

The  data  acquisition  program  was  simplified  when  high-resolution  reflection  was  dropped  from  consider¬ 
ation  as  an  effective  tool  for  imaging  acoustic  impedance  contrasts  at  depths  less  than  30  m  at  this  site. 
With  first  arrival  body  waves  and  surface  waves  left  as  the  only  seismic  energy  types  available  to  exploit 
at  this  site,  the  data  acquisition  requirements  were  reduced  enough  to  permit  recording  a  single  data  set 
optimized  for  both  analyses. 

At  most  sites  where  MASW  has  been  effectively  used,  low  natural  frequency  geophones  (4.5  Hz)  have 
provided  the  best  low  frequency  response  while  retaining  sufficient  high  frequency  components  of  the 
surface  wave  energy  to  produce  a  broadband  signal.  At  the  SWTR  site,  triple  10  Hz  geophones  provided 
the  best  response  for  the  depth  range  of  interest.  Additionally,  these  geophones  responded  well  enough  to 
first  arrival  body  wave  energy  for  automatic  first-break  picking  routines  to  operate  efficiently.  Optimal 
source  configurations  for  use  in  refraction  tomography  are  for  the  most  part  similar  to  those  employed  in 
surface  wave  recording.  It  is  generally  preferred  to  have  sources  that  are  high  energy  producers  with  a 
broadband  signal  centered  toward  the  lower  portion  of  the  conventional  seismic  spectra.  The  RAWD 
provided  a  broad  bandwidth  and  high-energy  pulse  that  produced  a  very  repeatable  waveform. 

To  insure  that  sufficient  offsets  were  recorded  for  refraction  tomography  and  that  close  trace  spacing  was 
maintained  for  MASW  analysis,  a  fixed  spread  of  240  receiver  stations  spaced  at  1 .2  m  was  deployed  for 
each  300  m  profile  (Figure  4).  A  240-channel  Geometries  StrataView  seismograph  was  used  to  record 
these  fixed-spread  data  using  a  source  spacing  of  4.8  m  inline  and  1.2  m  offline.  Each  source  location 
was  conditioned  with  an  initial  unrecorded  impact  to  seat  the  plate,  followed  by  three  recorded  impacts. 
Since  the  optimum  offset  window  for  surface  wave  sampling  of  the  upper  30  m  at  this  site  is  between  2 
and  25  m,  the  240  channel  spread  can  easily  be  split  into  appropriate  uniform  gathers.  Quality  control 
was  maintained  through  inspection  of  each  shot  gather  for  signal-to-noise  and  first  break  quality. 
Unacceptable  shots  were  deleted  and  re-recorded. 

High-resolution  elevation  data  were  acquired  with  accuracy  better  than  8  cm  (X,  Y,  Z).  This  level  of 
spatial  accuracy  was  critical  during  the  assignment  of  values  to  cells  within  the  study  volume.  If  future 
vehicle  tests  are  performed  at  SWTR,  these  spatial  data  will  provide  the  necessary  link  between  the  earth 
properties  and  sensor  deployments. 
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Figure  4.  Raw  shot  gather  with  first  arrival  picks  from  current  and  previous  shot  superimposed  on  section.  Low 
velocity  noise  arrivals  at  times  prior  to  the  first  arrivals  are  from  the  engine  of  the  weight  drop  source. 


7.  PROCESSING  OF  2V2-D  GRID 


7.1  Vs  Using  MASW 

Each  240-trace  shot  gather  was  split  into  a  left  and  right  spread.  The  optimum  60  traces  were  retained  to 
create  consistent  offset-distribution  shot  gathers  necessary  for  best  results  during  MASW  processing. 
These  60  trace  gathers  were  analyzed  with  SurfSeis  (a  proprietary  software  package  from  the  Kansas 
Geological  Survey  that  facilitates  use  of  MASW  for  continuous  profiling).  Each  shot  gather  generated 
one  dispersion  curve,  which  was  assigned  a  surface  location  corresponding  to  the  middle  point  of  the 
spread  being  analyzed  (Figure  5).  Care  was  taken  to  insure  that  the  spectral  properties  of  the  t-x  data 
(shot  gathers)  were  consistent  with  the  maximum  and  minimum  f-vc  values  (vc  is  the  phase  velocity  of 
surface  waves)  contained  in  the  dispersion  curve. 

SWTR  data  are  an  exception  to  the  general  “rule  of  thumb”  that  suggests  most  Rayleigh  wave  energy  is 
contained  in  the  fundamental  mode.  Instead,  higher  mode  Rayleigh  wave  energy  dominates  the  surface 
wave  frequency  band  on  most  shot  records.  A  variety  of  unique  processing  steps  was  undertaken  to 
eliminate  higher  mode  as  well  as  body  wave  energy.  A  two-phase  dispersion  curve  analysis  routine, 
requiring  extraction  of  the  fundamental  mode  portion  of  the  dispersion  curve  and  then  supplementing  this 
with  higher  frequency  portions  of  the  fundamental  mode  that  remain  after  time  domain  muting,  provides 
two  dispersion  curves  with  minimal  overlap  in  the  frequency  domain.  These  two  curves  were  digitally 
sutured  with  each  combined  dispersion  curve  individually  inverted  into  an  x-vs-z  trace.  In  general,  the 
combined  dispersion  curves  possess  a  useable  bandwidth  from  about  9  Hz  to  65  Hz.  Gathering  all  x-vs-z 
traces  into  shot-station  sequential  order  creates  a  2-D  grid  of  the  S-wave  velocity  field. 

To  increase  the  signal-to-noise  ratio  and  to  improve  convergence  during  inversion,  a  new  technique  called 
Filtering  of  Dispersive  Seismic  Event  (FDSE)  was  used  to  suppress  higher  modes  on  some  shot  records. 
This  method  removes  higher  modes  through  filtering  in  the  frequency  domain  and  avoids  the  detrimental 
artifacts  that  are  observed  in  dispersion  analysis  if  higher  modes  are  removed  using  time  domain  muting. 
Because  SurfSeis  was  designed  to  invert  fundamental  mode  Rayleigh  wave  energy  only,  inclusion  of 
higher  mode  or  body  wave  energy  can  inhibit  convergence. 
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Figure  5.  Data  example  from  key  output  portions  of 
the  acquisition  and  processing  flow,  a)  Raw  shot 
gather  edited  for  only  appropriate  trace  offsets  and 
propagation  directions,  b)  velocity  model  derived 
from  MASW  and  JASR  methods,  c)  Rayleigh  wave 
attenuation  coefficient  as  a  function  of  frequency  for 
these  data,  and  d)  Qp  and  Qs  inverted  from  attenua¬ 
tion,  frequency,  and  velocity  information. 
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7.2  Vp  Using  JASR 

Considering  the  resolution  requirements  and  redundancy  in  rays  penetrating  each  3mx3mx3m  cell, 
first  arrivals  were  picked  for  all  240  traces  on  every  fifth  shot  gather.  As  a  result,  1 7  shot  stations  pro¬ 
duced  the  Vp  cross-section  for  each  line.  Supplemental  analyses  demonstrated  that  only  minor  improve¬ 
ments  were  observed  (<  1 0%)  when  the  number  of  shots  processed  per  line  was  increased  to  63 .  As 
expected,  these  improvements  were  due  mainly  to  greater  ray  coverage.  However,  dropping  the  number 
of  shots  by  as  few  as  two  (down  to  15)  produced  noticeable  declines  in  data  uniformity  and  increased  the 
artifacts  associated  with  undersampling.  Ray  tracing  clearly  demonstrated  the  effect  of  the  sampling 
distribution  within  the  volume  in  general,  and  the  effect  that  each  cell  had  on  the  final  velocity  profile 
(Figure  6). 


Two-dimensional  Vs  cross-sections  obtained  from  MASW  analysis  were  used  to  generate  an  initial  inver¬ 
sion  model  for  the  tomographic  inversion  to  Vp  (Ivanov,  et  al.,  2000).  The  initial  model  was  optimized  by 
iterating  the  estimate  of  Poisson’s  ratio  until  the  first  arrivals  predicted  from  modeling  correlated  with  the 
actual  shot  records.  Several  inversion  runs  were  made  using  the  initial  model  to  converge  on  conditioning 
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Figure  6.  Ray  tracing  from  tomographic  analysis  of  first  arrival  data.  High  concentrations  of  rays  penetrate  all 
layers  of  the  model.  Especially,  high  densities  can  be  observed  through  layer  5  of  this  model. 


parameters  appropriate  for  this  data  set.  Fine-tuning  of  the  initial  model  is  optimized  when  best-fit  con¬ 
ditioning  parameters  were  used  during  preliminary  analyses. 

By  analyzing  the  correlations  between  the  model  and  observed  data  the  final  inversion  process  could  be 
used  for  quality  control  of  the  first  arrival  picking  routines.  In  some  instances,  secondary  first  arrival 
analysis  was  necessary  to  converge  on  a  “good”  solution.  Additional  quality  control  was  achieved  by 
verifying  that  the  2-D  Vp/Vs  data  were  reasonable. 

7.3  Q 

Estimates  of  Q  presented  in  this  paper  are  enhanced  developments  based  on  initial  work  described  by  Xia, 
et  al.  (2001)  and  unpublished  extensions  of  that  work  by  Xia  et  al.  Processing  these  data  to  obtain  Q 
involved  two  major  steps  for  the  first  two  methods  (A,  from  velocity  and  ar(f)  and  single  value  average  of 
near-surface  layer  based  on  A  approximation  with  calculation  of  Q  based  on  Eq.  A3).  The  third  method, 
inversion  of  Q  from  ar(f),  is  based  on  Eqs.  A4  through  A7  (Appendix  A).  Consistent  with  other  estima¬ 
tions  of  Q  in  the  shallow  subsurface,  Q  decreased  with  depth.  Q  has  been  demonstrated  to  have  a  wide 
range  of  responses  to  increases  in  depth  for  various  geologic  settings  (Xia,  et  al.,  2001).  Intuitively,  in¬ 
creasing  depth  should  increase  the  quality  factor  which  then  relates  to  decreases  in  attenuation.  However, 
regardless  of  the  method  used  to  calculate  Q  it  consistently  decreased  with  depth  in  the  upper  30  m  at  this 
site.  Decreases  in  Q  estimates  observed  when  using  approximation  methods  were  the  basis  for  develop¬ 
ment  of  the  three-method  approach  to  defining  Q  as  a  function  of  depth  necessary  for  model  simulations. 
More  extensive  and  detailed  model  simulations  will  eventually  provide  the  best  verification  for  these 
observations. 


8.  GEOSTATISTICAL  PROCESSING  OF  Vp,  Vs,  Qp,  AND  Qs 

Conventional  statistics  assumes  that  observations  used  in  estimation  are  independent;  that  is,  the  value  of 
one  observation  has  no  influence  on  the  value  of  a  subsequent  observation.  This  assumption  is  not  valid 
for  many  natural  phenomena  such  as  geophysical  properties,  because  a  value  at  one  location  obviously  is 
related  to  values  at  nearby  locations  and  less  closely  related  to  values  at  more  distant  locations.  Geo¬ 
statistics  is  a  special  branch  of  applied  statistics  that  takes  advantage  of  the  spatial  dependence  between 
samples  distributed  in  space  to  produce  optimal  estimates  of  a  property  at  unsampled  locations,  and  to 
assess  the  error  associated  with  these  estimates.  Although  the  geostatistical  literature  is  vast,  there  have 
been  relatively  few  applications  in  geophysics.  An  introduction  to  geostatistics  is  given  by  Davis  (2002) 
and  a  thorough  treatment  is  provided  by  Olea  (1999)  and,  from  a  statistical  viewpoint,  by  Cressie  (1993). 
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Information  on  spatial  continuity  is  derived  from  a 
model  of  the  semivariogram,  a  function  that  re¬ 
lates  distance  between  observations  to  the  mean 
difference  between  observations.  The  model  is 
fitted  by  weighted  least  squares  or  similar  proce¬ 
dures  to  empirical  values  of  semivariance  that  are 
calculated  from  the  observations.  For  a  property 
that  is  spatially  stationary,  the  semivariance  is 
inversely  proportional  to  the  autocorrelation  be¬ 
tween  observations.  A  theoretical  semivariogram 
(model)  is  used  rather  than  the  observed  empirical 
semivariogram  because  the  semivariances  must  be 
evaluated  for  any  distance,  not  just  the  discrete 
distances  between  observations  (Figure  7). 

The  estimator  used  in  geostatistics  is  kriging,  a 
family  of  generalized  linear  regression  techniques 
in  which  the  value  of  a  property  at  an  unsampled 
location  is  estimated  from  measured  values  at  neighboring  locations.  Kriging  estimates  require  prior 
knowledge  in  the  form  of  a  model  of  the  semivariogram  or  the  spatial  covariance.  Kriging  differs  from 
classical  linear  regression  in  that  it  does  not  assume  that  variants  are  independent,  nor  does  it  assume  that 
observations  are  a  random  sample. 

Kriging  is  used  to  make  regular  two-  or  three-dimensional  grids  of  estimates  of  geophysical  properties, 
which  constitute  the  required  earth  model.  Unlike  conventional  gridding  algorithms,  kriging  produces 
grids  that  have  statistically  optimal  properties.  Kriging  is  an  exact  interpolator;  that  is,  values  estimated 
by  kriging  will  be  exactly  the  same  as  the  observations  at  the  same  locations.  Kriging  estimates  are 
unbiased,  so  the  expected  value  of  the  estimates  is  the  same  as  the  expected  value  of  the  observations. 
Perhaps  most  importantly,  the  error  variances  of  kriging  estimates  are  the  minimum  possible  of  any  linear 
estimation  method,  and  the  error  variances  can  be  estimated  at  every  location  where  a  kriging  estimate  is 
made.  This  provides  a  way  of  expressing  the  uncertainty  of  a  property.  The  kriging  estimates,  in  combina¬ 
tion  with  the  estimation  variances,  can  be  interpreted  as  “lowest  likely”  values,  “most  likely”  values,  and 
“highest  likely”  values  of  the  geophysical  property  at  every  location  within  the  model  volume. 

In  general,  a  kriging  estimate  at  location  0  is  given  by  z0  =  Y'W  1 15  and  the  kriging  estimation  error 
variance  is<r0  =  B'W1  B,  where  Y  is  a  vector  of  the  observations  used  in  the  estimate,  W  is  a  matrix  of 
semivariances  between  these  observations,  and  B  is  a  vector  of  semivariances  corresponding  to  the  dis¬ 
tances  between  the  observations  and  the  location  being  where  the  estimate  is  to  be  made.  The  appropriate 
form  of  estimation  used  in  this  project  is  block  kriging,  the  general  name  for  a  collection  of  kriging  pro¬ 
cedures  that  include  a  change  from  point  observations  to  estimates  that  represent  the  average  of  volumes. 
Olea  (1999)  provides  a  succinct  development  of  block  kriging,  and  a  less  formal  discussion  is  given  by 
Isaaks  and  Srivastava  (1989).  In  block  kriging,  the  right-hand  vector  B  represents  the  averages  of  semi¬ 
variances  between  the  point  observations  x;-  and  all  possible  points  within  A,  a  volume  of  interest. 
(Actually,  the  elements  of  the  right-hand  vector  B  represent  the  semivariances  between  the  observations 
and  A,  integrated  over  the  volume  of  A.  As  a  practical  matter,  the  integration  must  be  approximated  by 
averages  of  semivariances  between  the  observations  and  point  locations  within  A.)  Typically,  determining 
the  kriging  estimate  and  estimation  variance  for  one  block  requires  solving  a  set  of  32  simultaneous  equa¬ 
tions.  Each  earth  model  based  on  the  first  15  lines  at  the  SWTR  site  contains  30,468,381  blocks  or  cells. 
There  are  two  parts  to  each  model,  one  containing  an  estimate  of  a  geophysical  property  such  as  Vp  in 
each  cell,  and  the  other  containing  the  estimation  variance  for  each  cell.  Models  have  been  prepared  for 
eight  geophysical  properties,  plus  a  2-D  model  of  surface  topography. 
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Figure  7.  Semivariogram  of  Vs  computed  along  seismic 
lines  in  the  SWTR  site,  YPG.  Open  circles  are  calcu¬ 
lated  semivariances;  fitted  line  is  a  best-fit  model 
(nested  Gaussian  model  with  a  range  of  60.7  m). 


To  produce  the  earth  model  of  a  geophysical  property  such  as  Vs 
(Figure  8),  we  first  compute  directional  semivariograms  and  examine 
them  for  nonstationarity.  If  nonstationarity  is  detected,  semivario- 
gram  models  are  computed  in  a  direction  perpendicular  to  the  drift 
and  used  for  block  kriging.  This  insures  that  a  conservative  set  of 
observations  will  be  used  in  calculating  the  kriging  estimates.  A 
series  of  experimental  models  are  fitted  to  the  experimental  semi¬ 
variances  and  the  model  producing  the  smallest  estimation  error  is 
chosen.  This  semivariogram  model  is  then  used  to  produce  block 
kriging  estimates  for  the  geophysical  property.  A  cross-validation  is 
performed  to  verify  that  the  block  kriging  estimates  are  valid  and  that 
the  observed  estimation  errors  are  within  the  theoretical  limits.  The 
resulting  three-dimensional  solid  model  can  be  displayed  as  cross- 
sections  along  any  desired  row  or  column,  and  as  maps  of  any 
desired  layer. 

The  results  of  geo  statistical  analyses  can  be  provided  in  the  form  of 
conditional  stochastic  realizations  (Deutsch  and  Joumel,  1998).  Each 
realization  is  a  three-dimensional  array  whose  values  conform 
exactly  to  the  observed  data  at  every  control  point  along  the  seismic 
lines.  At  all  interpolated  locations,  the  values  will  follow  the  same 
statistical  distributions  as  the  observed  data.  Each  realization  is  a 
“possible  scenario”  which  has  characteristics  that  are  statistically 
identical  to  those  actually  observed.  These  stochastic  alternatives  can 
be  used  to  evaluate  the  simulation  model  to  determine  the  model’s 
sensitivity  to  variation  in  geological  and  geophysical  properties. 
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Figure  8.  Shaded  contour  map  of  Vs 
at  5  m  depth,  SWTR  site,  YPG. 

Map  coordinates  are  in  meters  from 
an  origin  at  the  southwest  comer; 
contour  interval  is  50  m/sec. 


9.  SIMULATIONS  USING  3-D  VOLUMES  OF  Vp,  V„  Qp,  AND  Q 


Simulations  of  wave  propagation  through  a  KGS  Yuma  earth  model  were  performed  using  a 
order  finite  difference  staggered  grid  scheme  (Ptop)  and  directly  compared  to  field  data.  The 
variations  in  geophysical  parameters  and  topography  in  three  dimensions  make  seismic  wave 
a  highly  complex  phenomenon.  Successful 
synthetic  comparisons  with  field  data  give 
both  validation  of  the  Ptop  code  and  demon¬ 
strate  the  consistency  of  geologic  character¬ 
ization. 

The  preliminary  3-D  KGS  Yuma  model  had 
a  dimension  of  401  x  151  x51  nodes  with 
0.5  meter  spacing.  The  parameters  provided 
by  KGS  were  Vs,  Vp,  Qs,  Qp,  rho,  and  eleva¬ 
tion.  A  three-dimensional  view  of  Vs  is 
plotted  in  Figure  9.  The  Ptop  finite  differ¬ 
ence  code  allows  a  non-planar  surface  (Hest- 
holm,  1 999)  and  therefore  the  elevation  data 
(Figure  1 0)  was  incorporated  into  the  model. 

The  anelastic  (Q)  implementation  of  Ptop 
required  relaxation  time  constants  to  be 
calculated.  It  was  necessary  to  reduce  the 
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Figure  9.  KGS  3-D  S-wave  velocity  used  in  finite  difference 
calculations. 


computation  for  constants  from  40 1  x 
151  x51  xn,  where  n  is  the  number  of 
relaxation  mechanisms,  to  a  more  reason¬ 
able  number.  The  data  set  was  simplified 
by  rounding  velocities  to  the  nearest  tens 
of  m/s  and  binning  accordingly.  Model 
parameters  were  grouped  and  assigned  to  a 
representative  unique  velocity  block,  in 
this  case  91  unique  S-wave  velocities.  This 
reduced  computation  for  Q  relaxation 
parameters  to  91  x  n.  Computations  for 
both  the  elastic  and  anelastic  models  were 
run. 

The  time  domain  traces  were  normalized  on  the  maximum  amplitude  of  each  trace  to  emphasize  wave¬ 
form  characteristics,  arrival  times,  and  illustrate  the  relative  size  of  the  fundamental  mode,  the  particular 
phase  of  interest.  Time  domain  plots  of  the  normalized  velocity  amplitudes  for  the  anelastic  case  and  field 
data  are  shown  in  Figure  1 1.  The  synthetics  predict  accurate  arrival  times  for  both  the  P  and  surface 
waves  out  to  100  meters.  The  relative  amplitude  of  the  fundamental  mode  to  the  P  arrival  is  consistent 
with  the  data.  The  ringing  of  the  surface  wave  seen  in  the  field  data  is  not  reproduced  in  the  synthetics. 
This  may  be  due  to  fine  near  surface  layering  (less  than  2  meter)  where  the  model  is  poorly  constrained  or 
to  complexity  in  the  source  where  weight  bouncing  may  be  occurring. 

The  amplitude  decay  with  distance  emphasizes  the  agreement  of  signal  attenuation  with  distance  (Figure 
12).  Field  data  was  normalized  to  the  Ptop  Q  synthetic  at  the  16th  receiver  (approximately  20  m).  This 
was  approximately  the  range  where  the  clipping  of  field  data  stopped.  The  predicted  decay  rate  from  the 
synthetics  agrees  well  with  the  decay  of  the  field  data.  It  is  noted  that  the  Yuma  Q  model  used  here 
increases  with  depth. 
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Figure  10.  KGS  3-D  surface  elevation  used  in  finite  difference 
calculations. 
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Figure  11.  Record  section  of  Ptop  synthetics  (solid)  vs. 
field  data.  The  fundamental  features  of  P-wave  and 
surface  wave  arrival  times  and  relative  amplitude 
agree  well. 


Figure  12.  Amplitude  decay  vs.  distance.  Ptop  syn¬ 
thetics  (solid)  vs.  field  data.  Synthetic  calculated  with 
Q  is  normalized  with  field  data  at  approximately  20  m 
due  to  clipping.  Elastic  synthetic  decay  rate  is  plotted 
for  reference. 


Frequency  Power  Spectrum:  Ptop  Yuma  vs  Field  Data  @  91  m 
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Frequency  Power  Spectrum  @91  m  :  Yuma  Field  Data  vs.  Ptop  Synthetic 


Figure  13.  Frequency  spectrum  for  receiver  at  91.4 
meters.  Field  data  is  normalized  on  maximum  of  syn¬ 
thetic  computed  with  Q.  Spectral  content  and  slope 
match  for  30-110  Hz.  The  top  trace  is  the  elastic  syn¬ 
thetic  result  plotted  for  reference. 
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Figure  14.  Spectrum  of  field  data  and  synthetic  at  91  m. 
Synthetic  uses  Q  of  Yuma  model  and  includes 
instrument  response. 


Finally,  to  validate  the  frequency  content  of  the  synthetics,  the  spectra  from  receivers  at  91  meters  are 
plotted  in  Figures  13  and  14.  The  spectra  are  plotted  on  the  same  scale,  however  the  field  data  is  normal¬ 
ized  on  the  peak  amplitude  of  the  spectra  generated  by  the  Yuma  Q  model.  This  peak  frequency  corre¬ 
sponds  to  the  fundamental  mode.  The  implementation  of  Q  severely  attenuates  the  signal,  providing  peak 
frequency  match  as  well  as  the  slope  of  the  amplitude  from  30-110  Hz.  The  low  frequency  amplitude  that 
tapers  off  in  the  field  data  due  to  instrument  response  is  addressed  by  convolving  the  instrument  response 
with  the  synthetic  signal  (Figure  14).  The  synthetic  still  has  more  low  frequency  content  than  the  field 
data,  but  the  amplitude  does  drop  off  sharply  at  1 0  Hz.  The  slope  and  content  of  the  synthetic  spectra 
match  field  data  from  30  to  110  Hz.  There  is  also  somewhat  of  a  spectral  hole  at  55  Hz  that  is  matched  by 
the  synthetic  data.  The  maximum  spectral  power  at  30-45  Hz  in  both  the  field  data  and  the  synthetic  data 
corresponds  with  the  fundamental  mode  or  Rayleigh  wave. 

The  quality  of  agreement  is  extraordinary  and  rarely  achieved  in  practice.  This  is  due  to  the  fidelity  of  the 
geologic  model  and  the  successful  Yuma  site  characterization  performed  by  the  KGS/KU. 

10.  CONCLUSIONS 

A  high-resolution  geophysical  characterization  of  the  SWTR  site  at  Yuma  Proving  Ground  in  Arizona 
was  performed  by  the  KGS.  This  data  set  provides  valuable  insight  into  subsurface  features  that  may 
affect  future  tests  that  rely  on  seismic  wave  propagation  (such  as  unattended  ground  sensors).  By  putting 
constraints  on  subsurface  heterogeneities  and  surface  topography,  a  more  complete  understanding  of  local 
wave  propagation  may  be  achieved  leading  to  sensor  system  development  and  optimization. 

Synthetics  from  the  Yuma  model  with  Q  predict  the  amplitude,  dispersion  and  distribution  of  energy  in 
the  frequency  domain  observed  in  field  data.  The  agreement  is  a  preliminary  validation  of  finite  differ¬ 
ence  code  developed  at  ERDC-CRREL.  The  good  correlation  of  Ptop  synthetics  with  field  data  using 
simple  sources  bodes  well  for  future  applications  where  more  complex  sources  are  to  be  used.  The  ability 
to  predict  accurate  seismic  velocity  field  displacement  is  key  to  relying  on  synthetic  data  for  sensor 
system  development  and  optimization.  Reliable  synthetic  data  will  dramatically  reduce  system 
development  time  and  cost. 


11.  APPENDIX  A 


Consistent  with  all  three  was  the  calculation  of  the  attenuation  coefficient  as  a  function  of  frequency  at-.  The 
attenuation  coefficient  is  defined  by  the  following  expression: 


A  (  x  +  dx  )  =  A  (x)  e  "adx. 


(Al) 


Rayleigh  wave  amplitude  is  A,  a  is  the  Rayleigh  wave  attenuation  coefficient,  and  x  and  dx  are  geophone  location 
and  station  interval,  respectively.  After  the  FFT  with  respect  to  time,  the  expression  becomes: 


W(x  +  dx,f) 

x  +dx 

W(x,f) 

V  x 

(A2) 


Where  Otf  is  the  Rayleigh  wave  attenuation  coefficient  as  a  function  of  frequency  f  and  W  is  the  amplitude  at  a 
specific  frequency. 


After  calculation  of  at-,  the  processing  flow  diverges  for  each  of  the  three  methods.  The  first  two  methods  are  based 
on  the  assumption  of  a  homogeneous  medium.  The  third  method  is  based  on  a  layered  earth  model.  First,  Vik  from 
velocity  and  0Cf  requires  dispersion  analysis  to  establish  phase  velocity  of  the  Rayleigh  wave  relative  to  frequency. 
Calculation  of  Q(f)  then  follows  the  relationship: 

Q(f)=JtL  (A3) 

Vraf 

where  Vr  is  the  Rayleigh  wave  phase  velocity  as  a  function  of  frequency.  Moving  from  Q(f)  to  Q(z)  requires  Vr  and 
the  Vik  axiom  which  defines  the  surface  wave’s  depth  sensitivity  to  elastic  module.  Q(z)  in  this  case  is 
approximately  equivalent  to  Qs(z)  and  AQpfz). 

Second  is  the  assignment  of  a  single  Q  for  each  surface  shot  station.  This  method  follows  the  previous  flow  for  the 
Vik  approximation  method  except  Q(f)  is  averaged  for  all  frequencies  of  interest,  providing  a  single  Qs  and  Qp 
(based  on  the  Vi  relationship  to  Qs)  for  each  source  station. 


Finally,  by  inverting  at-,  Q  can  be  estimated.  The  relationship  between  Rayleigh  wave  attenuation  coefficients  and 
the  quality  factors  for  compressional  and  shear  waves  were  given  by  Anderson  et  al.  (1965)  as: 


(A4) 

where 

pi(f)=v11*'(t). 

’  dV„ 

(A5) 

s  (f ) = v  3C'W. 

,v  ’  S1  av. 

(A6) 

aR(f)  is  Rayleigh  wave  attenuation  coefficients  in  1/length,  and  f  is  frequency  in  FIz.  Qpi  and  Qsi  are  the  quality 
factors  for  compressional  and  shear  waves  of  the  ith  layer,  respectively.  Vp;  and  Vs;  are  the  P-wave  velocity  and  S- 
wave  velocity  of  the  ith  layer,  respectively.  Cr  is  Rayleigh  wave  phase  velocity,  n  is  the  number  of  layers  of  a 
layered  earth  model. 


Modifying  an  algorithms  discussed  in  Menke  (1984)  with  a  damping  factor,  it  becomes  possible  to  invert  for  Q 
using  the  following  statement  for  X;  >  0: 

AX  =  B  (A7) 

where  X  is  the  inverse  of  quality  factors  (model  vector  1/Q),  X;is  the  ith  component;  B  is  attenuation  coefficients 
(data  vector  a);  and  ,4  is  the  coefficient  matrix  determined  by  equation  (4). 

Figure  1.  Radial  walkaway  spread  deployed  to  investigate  the  antisotropy  of  geophysical  properties  at  the  SWTR 
site.  From  a  geostatistical  perspective  this  source/receiver  configuration  best  provided  estimates  of  geophysical 
directionality  in  the  test  area. 
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